In [1]:
%matplotlib inline
from pylab import *
from numpy import *
In [2]:
#dimensions of the square computational domain:
maxx = 10.
maxt = 1.
#difffusivity
D = 0.5
In [3]:
#discretization parameters
nx = 100 # number of unknown grid points in spatial direction
SC = 1.0 # stability criterion SC = 2*D*dt/dx^2, for FTCS should be SC <= 1
def diffusion_init(maxx, maxt, D, nx, SC, stride=1):
#choose time step according to CFL condition
dx = maxx/(nx+1)
dt = SC*dx**2/(4*D)
nt = int(maxt/dt)+1
#define array for storing the solution
U = zeros((int(nt/stride), nx+2, nx+2))
x = arange(nx+2)*dx
t = arange(nt)*dt
return U, dx, dt, x, t
In [4]:
def diffusion_solve(U, dx, dt, nx, nt, D, stride=1):
Ui = U[0,:,:]
for it in range(0,nt-1):
Ui[1:-1, 1:-1] = Ui[1:-1, 1:-1] + D*dt/dx**2*(Ui[2:, 1:-1] - 2*Ui[1:-1, 1:-1] + Ui[0:-2, 1:-1])\
+ D*dt/dx**2*(Ui[1:-1,2:] - 2*Ui[1:-1, 1:-1] + Ui[1:-1, 0:-2])
if it%stride == 0:
U[(it+1)/stride, :, :] = Ui
Let's try to propagate a square initial condition
In [5]:
U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=0.9)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
pcolormesh(U[0,:,:], rasterized=True, vmin=-1, vmax=1)
Out[5]:
In [7]:
diffusion_solve(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)
Out[7]:
The time step limit seems to be to strict, let's try increasing the time step twice
In [8]:
U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=1.1)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
In [9]:
diffusion_solve(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)
Out[9]:
In [ ]:
from scipy.sparse import dia_matrix, eye, kron
from scipy.sparse.linalg import factorized
def d2matrix(nelem):
elements = ones((3,nelem))
elements[1,:] *= -2
return dia_matrix((elements, [-1,0,1]), shape=(nelem,nelem)).tocsc()
def diffusion_solve_implicit(U, dx, dt, nx, nt, D):
alpha = D*dt/dx**2
mat2D = kron(d2matrix(nx), eye(nx)) + kron(eye(nx), d2matrix(nx))
M = eye(nx*nx)-mat2D*alpha
LU = factorized(M.tocsc())
for it in range(0,nt-1):
U[it+1,1:-1, 1:-1] = LU(U[it,1:-1, 1:-1].ravel()).reshape(nx, nx)
In [ ]:
U, dx, dt, x, t = diffusion_init(maxx, maxt, D, nx, SC=10)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
In [ ]:
diffusion_solve_implicit(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)
Again, Crank-Nicolson follows:
In [ ]:
def diffusion_solve_CN(U, dx, dt, nx, nt, D):
alpha2 = D*dt/dx**2*0.5
mat2D = kron(d2matrix(nx), eye(nx)) + kron(eye(nx), d2matrix(nx))
M1 = eye(nx*nx)-mat2D*alpha2
M2 = eye(nx*nx)+mat2D*alpha2
LU = factorized(M1.tocsc())
for it in range(0,nt-1):
U[it+1,1:-1, 1:-1] = LU(M2.dot(U[it,1:-1, 1:-1].ravel())).reshape(nx, nx)
In [ ]:
diffusion_solve_CN(U, dx, dt, nx, len(t), D)
pcolormesh(U[:, 30, :], rasterized=True, vmin=-1, vmax=1)
Neumann BC
In [ ]:
def diffusion_solve_CN_neumann_zero(U, dx, dt, nx, nt, D):
alpha = D*dt/dx**2
d2x = d2matrix(nx+2)
d2x[0, 1] += 1
d2x[-1, -2] += 1
mat2D = kron(d2x, eye(nx+2)) + kron(eye(nx+2), d2x)
M1 = eye((nx+2)**2)-mat2D*alpha/2
M2 = eye((nx+2)**2)+mat2D*alpha/2
LU = factorized(M1.tocsc())
for it in range(0,nt-1):
# U[it, boundary] -= alpha*dx*(NBC[it]+NBC[it+1])
U[it+1,:, :] = LU(M2.dot(U[it,:, :].ravel())).reshape(nx+2, nx+2)
In [ ]:
U, dx, dt, x, t = diffusion_init(maxx, maxt*50, D, nx, SC=50)
X, Y = meshgrid(x, x)
U[0, :, :] = 0.
U[0, (X<maxx*0.4) & (X>maxx*0.2) & (Y<maxx*0.4) & (Y>maxx*0.2) ] = 1.
In [ ]:
diffusion_solve_CN_neumann_zero(U, dx, dt, nx, len(t), D)
In [ ]:
pcolormesh(U[:, 30, :], rasterized=True, vmin=0, vmax=.1)
In [ ]:
# constants taken from http://en.wikipedia.org/wiki/File:Brusselator_space.gif
#initial field values
U0 = 1.
V0 = 1.
noise = 2.
#diffusion coefficients
DU = 0.2
DV = 0.02
#constant concentrations
A = 1
B = 3
def RHS(U, V):
fV = B*U - U**2*V
fU = A - fV - U
return fU, fV
maxt = 100.
maxx = 100.
U, dx, dt, x, t = diffusion_init(maxx, maxt, max((DU, DV)), nx, SC=.05)
V, dx, dt, x, t = diffusion_init(maxx, maxt, max((DU, DV)), nx, SC=.05)
X, Y = meshgrid(x, x)
U[0, :, :] = U0 + noise*randn(nx+2, nx+2)
V[0, :, :] = V0 + noise*randn(nx+2, nx+2)
#for i in range(10):
# U[0, randint(0, nx), randint(0, nx)] = 10
# V[0, randint(0, nx), randint(0, nx)] = 10
In [ ]:
def brusselator(U, V, dx, dt, nx, nt, DU, DV):
alphaU = DU*dt/dx**2
alphaV = DV*dt/dx**2
d2x = d2matrix(nx+2)
d2x[0, 1] += 1
d2x[-1, -2] += 1
mat2D = kron(d2x, eye(nx+2)) + kron(eye(nx+2), d2x)
M1U = eye((nx+2)**2)-mat2D*alphaU/2
M2U = eye((nx+2)**2)+mat2D*alphaU/2
M1V = eye((nx+2)**2)-mat2D*alphaV/2
M2V = eye((nx+2)**2)+mat2D*alphaV/2
LUU = factorized(M1U.tocsc())
LUV = factorized(M1V.tocsc())
for it in range(10):
U[0] = M2V.dot(U[it,:, :].ravel()).reshape(nx+2, nx+2)
for it in range(0,nt-1):
# U[it, boundary] -= alpha*dx*(NBC[it]+NBC[it+1])
# diffusion operator
U[it+1,:, :] = LUU(M2U.dot(U[it,:, :].ravel())).reshape(nx+2, nx+2)
V[it+1,:, :] = LUV(M2V.dot(V[it,:, :].ravel())).reshape(nx+2, nx+2)
fU, fV = RHS(U[it+1], V[it+1])
U[it+1] += fU*dt
V[it+1] += fV*dt
In [ ]:
brusselator(U, V, dx, dt, nx, len(t), DU, DV)
In [ ]:
figure(figsize=(12,12))
subplot(221)
pcolormesh(U[:, 30, :], rasterized=True, vmin=0, vmax=5)
subplot(223)
pcolormesh(V[:, 30, :], rasterized=True, vmin=0, vmax=5)
subplot(222)
pcolormesh(U[1400, :-2, :-2], rasterized=True, vmin=0, vmax=5)
subplot(224)
pcolormesh(V[1400, :-2, :-2], rasterized=True, vmin=0, vmax=5)
In [ ]: